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Figure 1. (a) Two atoms in state |+} collide and flip state in free space, (b) 
When these atoms are confined in a state-dependent optical lattice, their change 
of state must be accompanied by a tunneling event to a neighboring site. This 
leads to correlated hopping of atoms through the lattice. 

1. Introduction 

Ultracold neutral atoms arc a wonderful tool to study many-body physics and strong 
correlation effects. Since the achievement of Bosc-Einstcin condensation in alkali 
atoms [U El E], wc have been witnesses of two major breakthroughs. One is the 
cooling of fermionic atoms and the study of Cooper pairing and the BCS to BEC 
transition [4j [5l [6] . The other one is the implementation of lattice Hamiltonians using 
bosonic [71 [5] and fermionic atoms in off-resonance optical lattices. Contemporary and 
supported by these experimental achievements, a plethora of theoretical papers has 
consolidated ultracold atoms as an ideal system for quantum simulations. The goal 
is two-fold: AMO is now capable of implementing known Hamiltonians which could 
describe real systems in Condensed Matter Physics, such as Hubbard models [8] and 
spin Hamiltonians [21110]; but it is also possible to study new physical effects, such as 
the quantum Hall effect with bosons [TT] and lattice gauge theories [HI [13] . 

In this work we deepen and expand the ideas presented in [14j . where we 
introduced a novel mechanism for pairing based on transport-inducing collisions. As 
illustrated in Fig. [T] when atoms collide they can mutate their internal state. If 
the atoms are placed in a state-dependent optical lattice, whenever such a collision 
happens the pair of atoms must tunnel to a different site associated to their new state. 
For deep enough lattices, as in the Mott insulator experiments [7], this coordinated 
jump of pairs of particles can be the dominant process and the ensemble may become 
a superfluid of pairs. 

The focus of this work is to develop these ideas for all possible interaction 
asymmetries that can happen in experiments with bosonic atoms in two internal states, 
where the scattering lengths among different states can be different, {g^^ ^ g^) |15j . 
We want to understand the different types of correlated hopping that appear when 
we move beyond the limited type of interactions considered in Ref. |14j . We will 
better study the resulting phases and phase transitions, using different analytical and 
numerical tools that describe the many body states, and evolving from few particles to 
more realistic simulations. The main results will be to show that when one considers 
more general asymmetries, a new type of correlated hopping appears, but that both the 
previous |14| and the newly found two-body terms cooperate in creating a superfluid 
state with pair correlations. Indeed, this new state will also be shown to be more than 
just a condensate of pairs, based on analytical and numerical studies. 

Correlated hopping is not a new idea. It appears naturally in fermionic tight- 
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Figure 2. (a) Atoms in two internal states t and J, are trapped in an optical 
lattice and coupled by a Raman laser with Rabi frequency fi. (b) That setup is 
equivalent to two displaced superlattices for the dressed states |±) ~ |t) ± |J,). (c) 
When asymmetric contact interactions are considered and the hopping between 
superlattice cells neglected, the whole system behaves effectively as a ID array of 
alternating |+) and |— ) sites, with transitions of turmeling amplitude t between 
them. 

binding models, where it has been used to describe mixed valence solids |16j and, 
given that they are able to mimic the attractive interactions between electrons, also 
high- Tc superconductors [l71 [THl HH [201 HH [12 • In most of these works, the correlated 
hopping appears in the form UiO^-ak, indicating that the environment can influence the 
motion of a particle. This would seem substantially different from correlated motion 
of pairs of fermions. Nevertheless, even this more elaborate form of correlated hopping 
has been shown to lead to the formation of bound electron pairs [501 [10] and it has 
been put forward as a possible explanation for high Tc superconductivity |231 124] . 

This work is organized as follows. In Sec. [2| we introduce our model of correlated 
hopping ([1]), qualitatively discussing how it originates and what are the quantum 
phases that we expect it to develop. We present a possible implementation of this 
model which is based on optical superlattices and atoms with asymmetric interactions. 
Sec. [31 includes exact diagonalizations for a small number of atoms and sites. These 
calculations reveal the existence of insulating and coherent regimes, as well as pairing, 
and will be the basis for later analysis. In Sec. [4] we study the many-body physics 
of larger lattices with correlated hopping, using a variety of techniques: starting 
with insulating regime and following with the implementation of perturbation theory 
and the quantum rotor model. These methods suggest a number of possible phases, 
including a Mott insulator, a pair superfluid, a normal superfluid and a charge density 
wave state, and we estimate the parameters for which these phases appear. In Sec. [5] 
we develop two numerical methods to study our system, a Gutzwiller ansatz and an 
infinite Matrix Product State method. With these simulations we confirm the above 
mentioned phases and locate the quantum phase transitions, which are found to be 
of second order. Finally, in Sec. [5] we suggest some currently available experimental 
methods to detect and characterize these phases. 

2. Correlated hopping model 

We suggested in Ref. [H] that the combination of atomic collisions with optical 
superlattices can be used to induce correlated hopping. The basic idea is shown 
in Fig. [2b, where atoms are trapped in two orthogonal states called (+) and (— ). The 
interaction terms change the state of the atoms, forcing them to hop to a different 
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supcrlattice every time they collide. In this sense, interactions are responsible for 
transport. 

In this section we introduce the most general model of correlated hopping that 
can be produced by such means with two-states bosons. This model is presented in 
the following subsection, where we explain qualitatively the role of each Hamiltonian 
term. Later on, in Sec. 12.21 we establish the connection between the parameters of this 
model and the underlying atomic model. This is the foundation for the subsequent 
analytical and numerical studies. 

2.1. Lattice Hamiltonian 

In this work we study the ground state properties of a very general Hamiltonian that 
contains different kinds of correlated hopping. More precisely, the model will be 



Here, cj and Ci are bosonic operators that create and anihilate atoms according to 
the site numbering from Fig.Oa-c, and the colons : AiBj : denote normal ordering of 
operators Ai and Bj. 

Let us qualitatively explain the roles of the different terms in Eq. ([1]) . The first and 
second terms, U and V, are related to on-site and next-neighbor interactions. When 
these terms are dominant, we expect the atoms in the lattice form an insulator. Such 
a phase is characterized by atoms being completely localized to lattice sites, having 
well-defined occupation numbers, the absence of macroscopic coherence and a gapped 
energy spectrum. Whether this insulating state is itself dominated by strong on-site 
interactions U or by nearest-neighbor repulsion/ attraction V will decide whether it 
presents an uniform density, a Mott insulator (MI), or a periodic density pattern, a 
charge density wave (CDW), respectively. 

The third term is the key feature of our model. It describes the tunneling of pairs 
between neighboring lattice sites, with amplitude t. Given U, V, j = 0, we expect the 
atoms to travel along the lattice in pairs forming what we call a pair superfluid (PSF). 
These pairs will be completely delocalized, establishing long range coherence along 
the lattice. The observable (a^) would be the figure of merit describing this kind 
of delocalization, while a vanishing (a) indicates the absence of the single-particle 
correlations appearing in a normal superfluid. Furthermore, we expect this phase 
to have a critical velocity, similar to that of an atomic condensate, and the energy 
spectrum should be gapless. 

Unlike in Ref. [14j . when one considers the most general kind of atomic interaction, 
a second kind of correlated hopping appears, described by the last term in Eq. ([T|). 
Here, individual atoms will hop only if there is already a particle in the site they go to 
(cjcj(nj; — 1)) or leave at least a particle behind ((n^ — l)cjc|). One might be induced 
to think that this term is equivalent to single-particle hopping with a strength that 
depends on the average density, thus giving rise to a single-particle superfluid (SF) 
phase. However, this does not seem to be the case. We will show that the correlated 
hopping j generates a mixed phase which contains features of both the ordinary BEG 
and the PSF created by t. 





(1) 
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2.2. Relation to atomic parameters 

We now establish the relation between the model in Eq. ([1]) and the dynamics of atoms 
in an optical superlattice. The actual setup we have in mind is shown in Fig. ^a,- 
b and described in more detail in [Appendix A.l It consists on a three-dimensional 



lattice that is strongly confining along the Y and Z directions, creating isolated tubes. 
On top of this, we create an optical superlattice acting along the X direction [14]. This 
superlattice traps atoms in the dressed states |+) and |— ) , while the atomic interaction 
is diagonal in the basis of bare states | "f ) and | J,) . The interaction will be described by 
a contact potential and parameterized by some real constants gafs 



Hint = ^-f j €i^)4^^)M^)M^)dx (2) 



E 9al3 
2 

These interaction constants arc functions of the s-wave one-dimensional scattering 
lengths between different species g^p = 'i'Kh?a'^^p^ /m. In general, the interaction 
strengths among different atomic components are different from each other, a situation 
that can be enhanced with Feshbach resonances. We will use a parameterization 

9]i= 9a+ 91= gu, 9^1 = 90 + 92 , 911 = 90 ^92, (3) 
that makes the symmetries more explicit 

Hint = Y J ■ iPli^) + Pli^yf -(^^ + 9i J ■ P^ix)Pi{x) -.dx 

+ f j -.pAxf-Piixf-.dx (4) 

The total Hamiltonian combines the previous interaction with the kinetic energy and 
the trapping potential for one particle 

-^^' + V.{x) 



ips{x)dx, (5) 



^-Ey^i(-)L 2^ 

which is written in a different basis 

V'T(a;) = ^(V^++^'-), ^i(.T) = i=(^+-7^_). (6) 

Since the superlattice potential V±{x) is the dominant term, we may approximate the 
bosonic fields as linear combinations of the Wannicr modes in this superlattice and 
in the dressed state basis, a process detailed in [Appendix A[ Note that out of all the 
terms in the interaction Hamiltonian ([5]), only the first one is insensitive to the state 
of the atoms. This is important because the asymmetries gi and 527 when expressed in 
the dressed basis, produce terms that change the state of the atoms during a collision. 
Once we introduce the effective interaction constants in the lattice 

U,^ g, j \w{x)\'^dx fori ^0,1,2 (7) 

where w{x) is the single site Wannicr wavefunction, we arrive to the effective 

Hamiltonian in Eq. ^ with parameters V, t and j that relate to the microscopic 
model as follows 

2Uo + U, Ui U2 

c/ = ^^, ^ = -T' '=16 ^ = T- 

Unlike in the specific case of Ref. [Tl|, the most general situation contains not only 
two-body correlated hopping but also the terms proportional to j. 
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3. Preliminary analysis 

In this section we study the eigenstates of Hamiltonian ([T]) for systems that we can 
diagonahze exactly. The goals are to characterize the effect of the different interaction 
and hopping terms, as well as to understand the structure of the ground state 
wavefunction. Although we are limited to a small number of particles, the following 
examples provide enough evidence of the roles of correlated hopping, nearest neighbor 
repulsion and the utility of different correlators to characterize the states. 



3.1. A two-sites example 



Let us take the simplest interesting case: four particles in two sites. We write the 
Hamiltonian in the basis {|40) , |22) , |04) , |31) , |13)}, where the notation \n1n2) stands 
for ni particles in the first site and 712 in the second and we restrict to ni + ?i2 = 4, 

-4V6t -12j 
-6\/6j 


6V 
-12t 
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4/2 
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-4^/6^ 






-4V6t 
-6\/6j 





-4\/6t 



-12j 



\ 

-6^/6J 

-12j 
-12* 
W 



W. 



(9) 



Notice that in this particular case, U gives rise to a global energy shift and does 
not affect the different eigenstates. This is consistent with later studies where we 
will see that on-site interactions just add a global, density dependent contribution to 
the energy. To better understand the role of the remaining terms, we will consider 
separately three limiting cases, two of superfluid nature and an insulating one. 



Limit J 7^ 0, t = = 0, single-particle delocalization. In this case we take for 
simplicity = and diagonalize Eq. [51 finding the normalized ground state 

|Vo.*=o> = \ |40) + + \ 104) + \ 131) + \ |13) . (10) 

Note that this state is exactly a BEC of 4 particles spread over two sites 

|^BiJc(4)) = -^(^-^^ct^ 100). (11) 

This suggests that, at least in this small example, the correlated hopping proportional 
to j is equivalent to the single-particle hopping in the ordinary Bosc-Hubbard model, 
giving rise to the delocalization of individual particles. However, as it will become 
evident later on, for larger systems and more particles this interpretation is wrong. 



Limit j = 0, t ^ \V\, pair delocalization. In the presence of two particle hopping, 
the lowest energy state has the form 

|V'o,,=o) = C4o(t, V) |40) + C22(t, V) |22) + co4(t, V) |04) (12) 

with coefficients 

C22(i, 1^) OC -V^ + Vl2*2+-l/2, 

Cio{t, V) = co4(i, V) cx V6t. 



Correlated hopping of bosonic atoms induced by optical lattices 



7 



In particular, for dominant pair hopping t ^ |y| this is a state of defocahzed pairs 

1^) = I |40) + |22) + i |04) . (13) 

Observe that this wavefunction is not equivalent to what one would naively understand 
as a "pair condensate" from analogy with the single-particle case 

1^') ^ (E Ivac) - |40) + i |22) + 104) . (14) 

Instead the previous wavefunction is isomorphic to the BEG of two bosons 

|^Bi^c(2)) = i|20) + i=|ll) + i|02) (15) 

under the replacement of each boson with two atoms. It is also interesting to remark 
that V'B£;c(2) has larger pair-correlations than the state in Eq. ([M]). 



Limit \V\ — > oo, an insulator. Reusing the previous wavefunction (|12p and taking the 
limit of dominant nearest-neighbor interaction we obtain two possible states. For 
strong repulsion V +oo, states |40) and |04) are favored, forming a charge density 
wave (CDW) with partial fiUing \ipcDw) oc |40) -I- 104) . On the other hand, for strong 
nearest-neighbor attractions V —^ ~oo, the particles are evenly distributed forming a 
Mott insulator |22) . 



3. 2. Superfluidity of pairs 

We have seen that four particles in a two-sites lattice recreate the exact wavefunction 
of an ordinary BEG under the replacement of single bosons with pairs. We can test 
this idea for slightly bigger lattices, diagonalizing numerically the Hamiltonian which 
only contains the pair hopping term {t ^ 0, V, [/, j = 0). The resulting wavcfunctions 
are compared side by side with the BEG-like ansatz we mentioned. In the case of two 
particles we get indeed the expected result 

l^f-) = |V.r"')^X]«rivac), 

i=l 

whereas for 4 particles in 5 sites 

|i/,f -^ ) = ci (140000) + 104000) + 100400) -f |00040) 4- |00004) ) 

+ C2 (122000) + 102200) + |00220) + |00022) + |20002) ) 

+ C3( 120200) + 120020) + |02020) + |02002) + |00202) ), 

we find a disagreement between the ideal case of a BEG-like state with coefficients 
Ci = 1/5, C2 ~ c^ — \/2/5, and the exact diagonalization with ci ^ 0.2735, C2 ^ 0.3073 
and C3 ~ 0.1754. We observe that when compared to the ideal BEG, our paired state 
breaks the translational symmetry, revealing an effective attraction between different 
pairs, that favors their clustering. 

In Fig. [3] we plot the projection between these states, namely the solution of Eq. 
[l]with only t ^ and the ideal superfluid of pairs. In the nearby plot we also analyze 
two relevant correlators that will be also used later on in the manuscript, namely, a 
single-particle coherence 

^A = ^E('^lAa.) (16) 




Figure 3. (left) Fidelity between the ground state of the Hamiltonian and 
the ideal PSF for the case of 2 particles (line with circles) and 4 (no markers), 
(right) Particle fluctuation (solid), one-body correlator (dot-dash) and two-body 
correlator (dash) for the case of 2 particles (circle) and 4 (no markers) using the 
ground state of the Hamiltonian. 



and the pair correlator 

i 

As it is evident from the wavefunction and from the plots, there is no single-particle 
coherence or dclocalization because particles move in pairs. Hence, ^ Sao- The 
other correlator, C^, which we identify with the dclocalization of pairs is rather large 
and it only decreases with increasing the lattice size because the total pair density 
becomes smaller. 



4. Analytical methods 

We now study the many-body physics of our model for a much larger number of 
particles using exact analytical methods. We begin with the regime in which the 
interaction terms U and V dominate, obtaining the different insulator phases on the 
J, t = plane. Then, using perturbation theory, we compute the phase boundaries of 
these insulating regions for growing j and t. Finally, we study the properties of the 
ground state and its excitations in the superfluid phase, with j ~ and dominating 
t, proving indeed that this region describes a superfluid of pairs. 



No hopping limit: insulating phases 

To analyze the phase diagram it is convenient to work in the grand-canonical picture, 
in which the occupation is determined by the chemical potential fi. In this picture the 
ground state is determined by minimizing the free energy 

F = H -fiN, (18) 

where TV = Uk is the total number of particles, including both states |+) and |— ) . 
The free energy has a very simple form in the absence of tunneling 

4 

k '- 



F = 



{uk + Uk+if ■■ +Vnknk+i - fJ-Uk 



(19) 
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Figure 4. Phase diagram for U, V: different regions of stability and insulating 
phases for t, j = 0. 
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— (Hk + Uk+i )- + VukUk+i ~ I ) [nk + Uk+i ) 



e{x) 



This function is defined over positive occupation numbers Uk € {0, 1,2,.. .}. A discrete 
minimization will determine the different insulating phases and the regions where the 
system is stable against collapse. 

For a translational invariant system with periodic boundary conditions, all 
solutions can be characterized as a function of two integers a?* = (n, m), representing 
the occupations of even n2k = n and odd sites n2k+i = m. The optimization begins 
by noticing that the bond energy of two sites has a quadratic form 

C//4 {U + 2V)/i\ {U + 2fi)/4 

{U + 2V)/4: U/A J"" \ {U + 2fi)/4: 

= x^Ax-x^v, (20) 

where physical solutions are in the sector with rt, m > 0. For these occupation numbers 
to remain bounded, the bond energy e[x) has to increase as n, m or both grow. This 
gives us two conditions that need to be fulfilled to prevent collapse. If these conditions 
are not met, the ground state will be an accumulation of all atoms in the same site. In 
that case, the large interaction energies and the many-body losses induced by the large 
densities will cause the breakdown of our model and quite possibly of the experimental 
setup. 

The first stability condition is found by studying e(x) along the boundaries 
of our domain {n,m > 0). Take for instance m = 0, this gives a total energy 
£b ~ {U /'i)n^ — [{U + 2/x)/4]n. For this function to have a local minimum at finite n, 
we must impose 

U>0. (21) 

The second condition comes from analyzing the interior of the domain. For this, we 
take the only eigenvector of A which lies in the region of positive occupation numbers, 
n = m = x/2. This line has an energy = [{U + V)/4:]x'^ — [{U + 2fi)/A]x and, to 
have again a finite local minimum, we require 

V > -U. (22) 
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Given that Eq. (PT|) and Eq. (|22p arc satisfied, the system is stable and we have 
two possibilities to attain the minimum energy: either at the boundaries, n = or 
m = 0, or right on the eigenvector of A. Inspecting Eb and e+ we conclude that a 
positive value of V will lead to the formation of charge density waves (CDW) of filled 
sites alternating with empty sites 

y > ^ 7i2fc = or nsfe+i = 0. (23) 

If < our energy functional will be convex and the minimum energy state will be 
a Mott insulator with n ~ m, when n + m is even, or a charge density wave with 
n = m± 1, when n + m is odd. The actual choice between these two insulating phases 
is obtained by computing the energy of both states 

e(2n+ 1) = U{2n+lf/A + Vn{n+l) - (?/ + 2/i)(2?7, + l)/4 (24) 

e(2n) =Ui2nf/4: + Vn^-iU + 2^i)2n/A. (25) 

Having e(2n+ 1) — e(2?i) = defines the value of fj. at which the state with 2n particles 
every two sites, a Mott with n particles, stops being the ground state and becomes 
more favorable to acquire an extra particle to form a CDW. The boundaries of these 
insulating phases for t, j = arc given by 

li{2n^2n + l) = {U + V)2n, (26) 

li{2n-l ^2n) = {U + V)2n-U. (27) 

Thus summing up, for /x(2n— 1 2n) < fJ, < /i(2n 2n+l) the optimal occupation is 
n particles per site, forming a Mott, while for ij,{2n 2n+l) < ji < /x(2n+l — > 2n+2) 
the occupation number is 2n + 1 particles spread over every two sites, having a CDW. 
The results of this section are summarized in Fig. 14.11 



^.2. Perturbation theory: insulator phase boundaries 

The previous calculation can be improved using perturbation theory for t, j U,V 
around the insulating phases, obtaining the phase boundaries around the insulators 
as t and j are increased. This is done applying standard perturbation theory up to 
second order on both variables [25j . using as unperturbed Hamiltonian the operator 
P9p and as perturbation the kinetic energy term 



W 



{- ^^c'+i + J {n. - l)cl(c,-i + c,;+i)] + H.c.} . (28) 



We start calculating analytically the ground state energies of the first four 
insulating phases according to ([T9)) . considering the perturbation W up to second 
order in t. For the CDW with Ui = 1 and n^+i = this energy is obviously zero 

E{L/2) ^ 0. (29) 

For the MI with one particle per site we have virtual processes of the correlated 
hopping j, as environment-assisted hopping starts being allowed in an uniformly filled 
lattice 

EiL) = iU + 2V)L/2-jj^L. (30) 

For the CDW with Ui = 2 and 71.;+ 1 = 1, we find some doubly occupied sites and 
contributions from the pair hopping t 

, . , . . fQt^ 24j' 32^2 \ 

EiL + L/2) = iSU + W)L/2 - (- + ^ + ^ j L. (31) 
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Figure 5. Phase diagrams for fi, j,t ^ and V = —0.05. (left) Varying t 
for J = 0. (right) Varying j for t = 0. In both cases the lowest region is a CDW 
with alternating and 1 particle occupation, followed upwards by a Mott of one 
particle per site, a CDW with 1 and 2 particles and the highest area a Mott of 
two particles. 



Finally, for the MI with two particles per site, a calculation detailed in 



E{2L) = {3U + AV)L 



-L. 



(32) 



U-2V 

Here L is the total number of sites and all results presented in this section are for the 
case V < 0. 

At each value of j, t, the boundary of an insulating phase with average density n 
is given by the degeneracy condition with a compressible state E{nL) = E{nL ± 1). 
Those points correspond to the chemical potential at which a hole, fj,h{nL) ~ 
E{nL) — E{nL—l), or a particle, iip{nL) = E{fiL + l) — E{nL), can be introduced We 
show here the lower and upper limits of the first four insulating regions, corresponding 
to the CDW with = 1, n.i^i = 

l^h{L/2) = 
fip{L/2) =U + 2V + 
the Mott with one particle per site 

HhiL) = U + 2V + f 
fip{L) ^E{L + l 

^2U + 2V - 8j - 
the CDW with n.; = 2, Ui+i = 1 

^ih{L + L/2) ^2U + 2V 

'4 



V U' 



4 
U 
E{L) 



2 
V 



16 



U -2V 



4f 
U 



2f 

V ' 

- 2UV 
2 96 



8{f - 6jt - 3t^) 
U-2V 



(33) 



(34) 



(35) 



24 



128 



fip{L + L/2) = E{L + L/2+l) 



U-6V U 
-E{L + L/2) 
108 54 



2V U + 2V 



(36) 



96 
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3U~2V U-2V U~W, 



(37) 



and the MI with two particles per site 
Aih(2L) = 3C/ + 4F + 8 



108f 
U 



2U^ 



54j2 24(35^2 - 3\/2ji 



U 



^p(2L) 



V 



32(57j2 



U -2V 
lOSf 



A8t^ 



6V 3U-2V 



48^2 



U U-6V 
12^/6j^ - 7t^) 



3U-2V 



(38) 



(39) 



V 3U~6V 
The corresponding boundaries are plotted in Fig. [5l For small hopping amplitude, 
they match the values that are found later on with the numerical methods. But even 
for larger values, this approximation anticipates that the lobes are significantly larger 
for pair hopping t than for the correlated hopping j. 

4-3. Phase model: analysis of the pair condensate 

So far we have studied the many-body physics around the limit of strong interactions. 
However, the main goal of this work is to understand the effect of correlated hopping 
and the creation of a pair superfluid. In absence of a mean field theory, but still 
in the limit of dominant two-body hopping U,V ^ t, we can use the number-phase 
representation, introduced in Ref. [27] for an ordinary BEC. Note, however, that the 
model in Ref. [27] cannot be directly applied here. Following that reference, one would 
assume a large number of particles per site, Ui > 1, and introduce the basis of phase 
states \<f>) 



{ri\$) = (27r)-^/2e™-f 



(40) 



Using these states, one would then develop approximate representations for the 



operators af 



.t2 



and Hi, and diagonalizc the resulting Hamiltonian in the limit of 



weak interactions. But after a few considerations one finds that the resulting phase 
model does not preserve an important symmetry of our system: if j = particles can 
only move in pairs and the parity of each site, (—1)"', is a conserved quantity. 

To describe correlated hopping we must use a basis of states with fixed parity 



{2n + iy'\^), = (27r) 



- L /2 ^in-i 



iye{0,l} 



(41) 



which is v ~ for the ground state we are interested in. As mentioned before, we now 
have to find expressions for the different operators, af, aj^ and rii. We use the fact 
that our states will have a density close to the average value n and approximate the 
action of the operators over an arbitrary state as 

(42) 
(43) 

.t2 ^\ _ ^ /7irT~rT7i~roT„i</'. A (-44) 



$) ^ ^n{n-l)e-''l'' \$), 
$) = v/(" + l)(^ + 2)e*^' |<?) 
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Introducing the constant = n{n—l){n + l){n + 2) our Hamiltonian becomes similar 
to the quantum rotor model |27| 

L 

H^Y. + (2t^ + ^"^^^A^^ + 2p'*cos (0, - 0,+i)] 

1=1 

- (C/ + l/)Lfi2. (45) 

For small U and V, the ground state of this model is concentrated around — 0^+1 = 0. 
Expanding the Hamiltonian up to second order in the phase fluctuations around this 
equilibrium point, we obtain a model of coupled harmonic oscillators. This new 
problem can be diagonalized using normal modes that are characterized by a quasi- 
momentum k = 27rn/L, n € [-{L - l)/2, (L - l)/2]. The result is 

H = y^hLOk (hlhk + ^+Ea (46) 



with normal frequencies 



Uk = ipVmt \sm{k/2)\ ^1 + (^1 + cos (fc), (47) 

and a global energy shift Eq = 4(C/ + V)N^ - 4:Ln^{U + V)- 2pHL. 

It is evident from Eq. (j47|) that our derivation is only self consistent for negative 
values of V. Otherwise, when V > some of the frequencies become imaginary, 
signaling the existence of an unbounded spectrum of modes with > 7r/4 and that 
our ansatz becomes a bad approximation of the ground state. This strictly means 
that our choice (jji = only applies in the case of attractive nearest neighbor 

interactions, —U < y < 0, as we know that this interaction cannot destabilize a 
translational invariant solution such as the uniform Mott insulator. However, it does 
not mean by itself that the whole system becomes unstable for V > — indeed, we 
will show numerically that it remains essentially in a similar phase for all values of V, 
but in the case of > the insulating phases are stable until values of the hopping 
slightly higher as in the < case. 

If we focus on the regime of validity, we will find that the spectrum is very 
similar to that of a condensate. At small momenta the dispersion relation becomes 
linear, LUk (x Vgk, with sound velocity Vg = Ap^/2Ut/h^ while at larger energies the 
spectrum becomes quadratic, corresponding to "free" excitations with some mass. 
This a consequence of the similarity between our approximate model for the pairs ()45p 
and the phase model for a one-dimensional condensate. However, we can go a step 
further and conclude that the similarity extends also to the wavefunctions themselves, 
so that the state of a pair supcrfluid can be obtained from that of an ordinary BEG 
by the transformation n s- 2n. This is indeed consistent with what we obtained for 
the diagonalization of a two-particle state in the limit j, U,V = [See Eq. (fT5|l ] . 



5. Numerical methods 



The previous sections draw a rather complete picture of the possible ground states in 
our model. In the limit of strong interactions we find both uniform insulators and a 
breakdown of translational invariance forming a CDW, while for dominant hopping 
we expect both single-particle superfluidity and a new phase, a pair superfluid. We 
now confirm these predictions using two different many-body variational methods. 
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5.1. Gutzwiller phase diagram 

The first method that we use is a variational estimate of ground state properties based 
on a product state [15] 

|^GW>=nE/"^7=T"I"|0). (48) 

Minimizing the expectation value of the free energy F = H — i^iN with respect to the 
variables /„, under the constrain of fixed norm |/nP = 1, we will obtain the phase 
diagram in the phase space of interactions and chemical potential {U, V, j, t, /i). 

In our study we have made several simplifications. First of all, we assumed 
period-two translational invariance in the wavefunction, using only two different sets 
of variational parameters, fn^'^^^ — and f^^^ = f^. In our experience, this is 
enough to reproduce effects such as the CDW. Next, since J7 > is required for the 
stability of the system, we have taken [/ = 1 as unit of energy. The limit C/ = is 
approximated by the limits j, i 3> 1 in our plots. Finally, in order to determine the 
roles of j and we have studied the cases j = and t = separately. The results are 
shown in Fig. [6] and Fig. [7] for V <Q and ^ > 0, respectively. 

The first interesting feature is that, as predicted by perturbation theory, we 
have large lobes both with integer 1,2,... and with fractional 1/0, 2/1, . . . occupation 
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Figure 7. Phase diagram witli a GW ansatz for U = 1,V = 0.5. Wc focus 
on (left) j = and (right) t = 0, separately. Upper plots show average site 
occupation (n) (grayscale) and number variance An (contour). Lower plots show 
single-particle coherence (a) (grayscale) together with Aa^ (contour). 



numbers, forming uniform Mott insulators and CDW, respectively. The insulators 
are characterized by having a well defined number of particles per site, and thus no 
number fluctuations An^ = ("'^) ^ ("-) ^ = 0. While the size of the lobes does not depend 
dramatically on the sign of V, these are significantly larger for the pair hopping t than 
for the correlated hopping j, as already seen with perturbation theory. 

The boundary of the insulating areas marks a second order phase transition 
to a superfluid regime, where we find number fluctuations An ^ 0. In order to 
characterize these gapless phases we have computed the order parameter of a single- 
particle condensate (a), and two quantities that we use to detect pairing. The first 
one is a two-particle correlation that generalizes the order parameter of a BEC to the 
case of a pair-BEC (a^). The second quantity Aa^ = |(a^) — (a)^| is used to correct 
the previous value eliminating the contribution that may come from a single-particle 
condensate coexisting with the pair-BEC. 

When j = we always find that (a) ~ 0, even outside the insulating lobes. This 
marks the absence of a single-particle BEC, which is expected since we do not have 
single-particle hopping. On the other hand, we now find long range coherence of the 
pairs and thus (a^) 7^ all over the non-insulating area, which we identify with the 
pair-superfluid regime. 

The situation is slightly different for t = 0. The single-particle order parameter 
(a) no longer vanishes in the superfluid area, denoting the existence of single-particle 
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coherence, but at the same time we find that the two-particle correlations exceed 
the contribution from the single-particle superfluid as Aa^ 7^ 0, which we attribute 
to a coexistence of both a single-particle and a pair-superfluid, or a state with both 
features. 

This picture does not change substantially when V is positive or negative. The 
only differences are in the insulating regions, where the CDW is either due to the 
incommensurability of the particle number {V < 0) or really gives rise to the separation 
of particles alternating holes and filled sites {V > 0). However, in the superfluid regime 
we find no significant changes and in particular we see no breaking of the translational 
invariance or modulation of the coherent phase. 



5.2. Matrix Product States: long range pair correlations 

The previous numerical simulations are very simple and cannot fully capture the 
single particle and two-particle correlators. To complete and verify the full picture we 
have searched the ground states of the full Hamiltonian using the so called iTEBD 
algorithm, which uses an infinite Matrix Product State ansatz together with imaginary 
time evolution 

Roughly, this ansatz is based on an infinite contraction of tensors that 
approximates the wavefunction of a translational invariant system in the limit of 
infinite size. Adapting the ansatz to our problem we write it as 



n 



X , ^, A lti^'^llir Ivac> . (49) 

Here the r° and are matrices that depend on the state of the odd and even sites 
they represent, a dependence which is signaled by the 712^+1 and 7i2fc+2 hi the previous 
equation. These matrices are contracted with one-dimensional vectors of positive 
weights A^'° > 0, which are related to the coefficients of the Schmidt decomposition. 
This variational ansatz is known to work well for states with fast decaying correlations, 
but it also gives a good qualitative description of the critical phases. 

In order to optimize the iTEBD wavefunction we performed an approximate 
imaginary time evolution using a Trotter decomposition and local updates of the 
associated tensors, as described in Ref. [30]. Using the canonical forms for these 
tensors it is also straightforward to compute expectation values for different operators 
acting either on neighboring or separated sites. 

In Fig. [5] we plot the most relevant results for three cuts across the phase diagram, 
II = 0.5, 1.5 and 2.5, so that each line crosses both an insulating plateau and the 
superfluid region. We have used small tensor sizes from Z? = 16 up to 64, a value 
limited by the need of using large cutoffs for the site populations {n„iax =8). As 
shown in the figures, when j = the single-particle correlator is zero for distinct sites, 
and we are left only with two-particle correlations. In the MI case the pair correlations 
between neighboring sites decrease very quickly, while in the superfluid regime we see 
a critical behavior 

Cl « (50) 

with an exponent that varies between a = 0.5 and a ~ 0.6, depending on the 
simulation parameters. 
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Figure 8. Numerical simulations with the iTEBD algorithm. On the left we 
plot the results for j = 0, (7 = 1, V = —0.05 and on the right for V = 0.05, as a 
function of two-particle tunneling t. The upper row shows the density (dashed) 
and particle number fluctuations (solid), while the lower row shows (Aa)^ (solid) 
and two-particle coherence (dashed), which overlap indicating (a) = 0. For each 
set of operators we show three lines, for ^/f/ = 0.5, 1.5 and 2.5 corresponding to 
plain line, circle and star, respectively. 



6. Detection 

There arc many ways of differentiating the phases we have found, each one having 
its own degree of difhculty. The simplest and best estabhshed detection methods are 
related to the insulating regimes. These phases, which involve both MI and CDW, are 
characterized by having a well defined number of particles at each lattice site, the lack 
of coherence and an energy gap that separates the insulator from other excitations. 

The energy gap in these insulators may be probed either by static or spectroscopic 
means as it has been done in experiments [71 131j , determining that indeed the system 
is insulating. Second, the lack of coherence will translate into featureless time-of-flight 
images, having no interference fringes [7] at all. Even though there will be no fringes, 
the measured density will be affected by quantum noise. The analysis of the noise 
correlation will show peaks at certain momenta [3ll[33| that depend on the periodicity 
of the state, so that the number of peaks in the CDW phase will be twice those of the 
MI. In case of having access to the lattice sites, as in the experiments with electron 
microscopy |34| . or in future experiments with large aperture microscopic objectives 
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Figure 9. Correlator = vs. site separation, for {t, fi) = (0.1,1.5), 

(0.3, 0.5), (0.2, 1.5), (0.1, 2.5), from top to bottom. The line with t = 0.1, = 1.5 
(solid, circles) corresponds to a MI state and correlations decay to zero on the 
third site. The remaining lines have been rescaled and plotted in log-log scale in 
the inset, which shows that the two-particle correlator decays approximately as 
with a power a ~ 0.5 — 0.6. 



that can collect the fluorescence of individual lattice sites [35l [36] , the discrimination 
between the MI and the CDW should be even easier, since in one case we have a 
uniform density and in the other a periodic distribution of atoms. 

When the system enters a supcrfluid phase, it becomes a perfect "conductor" 
with a gapless excitation spectrum. The lack of an energy gap, should be evident in 
the spectroscopic experiments suggested before. However, we are not only interested 
in the superfluid nature, but rather in the fact that this quantum phase is strongly 
paired. More precisely, we have found that for j = the single-particle coherence is 
small or zero, and that the two-particle correlator decays slowly 



The first equation implies that the time of flight images will reveal no interference 
fringes and will exhibit noise correlations which will be similar to the MI. In order to 
probe and confirm the pairing of the particles, we suggest to use Raman photo- 
association to build molecules out of pairs of atoms [371 [35]. For an efficient conversion, 
it would be best to perform an adiabatic passage from the free atoms to the bound 
regime. As described in |39], we expect a mapping that goes from \2n) — > |n) , where 
2n is the number of bosons and n the number of molecules. More precisely, we expect 
the a^^ operator to be mapped into , so that the pair coherence of the original 
atoms translates into the equivalent of for the molecules. This order should reveal 
as an interference pattern in timc-of-flight images of the molecules. 

Finally, in cases with j ^ 0, we have found the coexistence of single-particle 
and two-particle coherences. This translates also into the coexistence of interference 
fringes with nonzero pair correlators. 




(51) 
(52) 
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7. Conclusions 

Summing up, we have suggested a family of experiments with cold atoms that would 
produce correlated hopping of bosons. The mechanism for the correlated hopping is 
an asymmetry in the contact interaction between atoms. This asymmetry is exploited 
by trapping the atoms in dressed states, a configuration that gives rise to transport 
induced by collisions. 

The main result of this paper is that there is a huge variety of interaction 
asymmetries that will give rise to long range pair correlations via interaction-induced 
transport. Formally, in the resulting effective models we recognize two dynamical 
behaviors. If we have a nonzero asymmetry in the interspecies interactions, 

t « 2.911 - .git - gii ^ 0, 

the Hamiltonian will exhibit pair hopping, while an asymmetry in the intra-species 
scattering lengths 

j « ffTT - 911 7^ 0, 

gives rise to correlated hopping. However, we have given enough evidence that both 
Hamiltonian terms give rise to a novel quantum phase which we call pair supcrfluid. 
This phase is characterized by a gapless spectrum with a finite sound speed, zero 
single-particle correlations and long range pair coherence. All quantum phases are 
connected by second order quantum phase transitions. These phases can be produced 
and identified using variations of current experiments [40U411I55] . The nonperturbative 
nature of the effect should help in that respect. 

Our ideas arc not restricted to one dimension. It is possible to engineer also a two- 
body hopping using two-dimensional lattice potentials. Again, the basic ingredients 
would be atoms with an asymmetric interaction and an optical lattice that traps two 
states, 1+) and |— ) with a relative displacement. Both in the one- and two-dimensional 
cases it is a valid question to ask whether the coupling between different trapped 
states, |±) , can excite also transitions to higher bands, processes that have not been 
considered in the paper. Our answer here is no. There arc only two sources of coupling 
to higher energy bands. One is the interaction, but we arc already assuming that the 
interaction energies are much smaller than the band separation. Following the notation 
from Ref. [8], we have the constraint that the interaction energy should be smaller 
than the energy separation to the first excited state in a well of the periodic potential, 
fi^U ^ hvn the same requisite as for ordinary Bose-Hubbard models [7]. The other 
source of coupling to higher bands would be single-particle hopping. However, unlike 
[5], here we arc assuming that these terms arc strongly suppressed compared with 
the interaction. In other words, realizing the models that we suggest in this paper, 
for realistic densities, n = 2, and simple potentials, imposes no further constraint in 
current experiments. 

Finally, let us remark that transport-inducing collisons may be implemented using 
other kinds of spin-dependent interactions. For instance, correlated hopping appears 
naturally in state-dependent lattices loaded with spinor atoms, because their inter- 
actions can change the hyperfine state of the atoms while preserving total angular 
momentum [42] . 
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Appendix A. Derivation of the model in superlattices 

As discussed before, the main idea behind atomic correlated hopping is to trap atoms 
whose interaction allows them to change their state. In this section we provide one 
possible implementation of this idea using state dependent superlattices that trap 
dressed states. 



Appendix A.l. Dressed states trapping 

Our starting point is the setup in Fig. [^Ji, which is itself taken from Ref. [13]. It 
consists on an optical lattice trapping atoms in states |t) and ||), together with a 
Raman coupling between these states. Mathematically, this configuration is described 
by the single-particle Hamiltonian 

Htr.p = Vo sm{kxf{\^) (T| + U)ai) 

+ f] sin(fcx)(|T)ai + U)(T|). (A.l) 
By moving to the basis of dressed states |±) = (||) ± ||)) , we find that the trapping 
is effectively equivalent to two superlattices with a relative displacement as in Fig. [2h 

i?trap = {Vosm{kxf + nsm{kx)) 1+) (+1 

+ (Vbsin(fcx)2 -nsm{kx)) |-) (-| . (A.2) 

Under appropriate circumstances, discussed in [43| . we will find that each supcrlattice 
site has a unique ground state, energetically well differentiated from the next excited 
state, and which consists of a symmetric wavefunction spanning both lattice wells. If 
this is the case and if all energy scales, such as the interaction and the hopping, are 
small compared to the separation between Bloch bands, we can expand the bosonic 
field operators describing the atoms in terms of these localized wavefunctions 

= ^C2. iy(x-2i0, (A.3) 

i 

i,. (x) = ^2^+1 - + 1)0, 

i 

where 21 = 27r//(; is the supcrlattice period, Cj are bosonic operators that, for j even 
(odd), annihilate an atom in state 1+) (|— )) in the j-th superlattice cell 

C2i = ;^('^2i,+ + a2i+i,+), (A. 4) 

1 , , 

C2i+1 — -^(a2i+l,- +a2i+2,-), 

and the localized wavefunctions Wix) are a superposition of the Wannier functions 
w(x'^ of the underlying lattice 

W{x - 2H) ^ ^ \w{x - 2H) + w{x - (2i + 1)/)] . (A.5) 
v2 
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Appendix A. 2. State- changing collisions 

We will now express the interaction @ in the basis of dressed states. We proceed 
using the change of variables in Eq. [5] to find the expression of the densities 

ptW = ^(p+ + p- + ^+^- + 'A1V'+), (A.6) 

Pl{x) = \{p++P~-^li^--ijU^+)- (A.7) 

The first obvious conclusion is that the total density is independent of the basis on 
which it is written, 

p(.t) = (.t) + {x) = p+ {x) +p-{x). ( A.8) 

Hence, the term of g^ is insensitive to the state of the atoms. On the other hand, the 
asymmetric terms arc not so simple. The gi interaction, which is proportional to the 
product of densities 

: P^Pl ■■ = \-[p++ P-f ■.-\- (VAV'- + : 

= i:p^+p2_:_l(v,t_V^+i7.e), (A.9) 

gives rise to a scattering that changes the state of interacting atoms from |— ) to |+) 
and viceversa, as in Fig. [TJi. The term of 92 has a lightly different effect, 

:pT(a;)'-Pi(x)2 : (A.IO) 

it gives rise to processes where one atom changes its state influenced by the surrounding 
environment. In the following subsections we will see what happens to the interaction 
terms (jA.8|) . (|A.9p and (jA.lOp . when the atoms are confined in a lattice. 

Appendix A. 3. Final model 

In this section we will put the previous results of this appendix together. We will 
take the tight-binding expansion of the field operators (|A.3p and use it together 
with Eqs. (|A.8|) . (jA.9[) and (jA.lOp to expand the interaction Hamiltonian (|4]). For 
convenience, we will rename the bosonic operators as 

C2k = ak+ and C2k+i = (A. 11) 

according to the position at which their Wannier functions are centered (see Fig. 
Along the derivation, one obtains many integrals of ground state wavefunctions 

Ck,„,= J \W{x- kl)\^\W{x-ml)\^dx. (A.12) 

We will only keep those integrals with a separation smaller than a superlattice period. 
Taking Eq. (|A.5[) . the expression for the superlattice localized states, one obtains 

Ck,k ^ f \Wix)\Ux^l f \wix)\Ux, (A.13) 



Ck,k±i = J \W{x)\^\W{x - Z)pdx ~ i y \w{x)\Ux, (A.14) 
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where wlx) are the Wannier wavefunctions of the underlying sublattice. Using these 
tools, the symmetric interaction term becomes 



90 
2 



dx : {p^{x) + Piix)y 

N/2 



^ ^ : n\f.C2k,2k + n2kj^iC2k+l.2k+l + 2n2kn2k+lC2k,2k+l 



go 

4 



N 



■ / dx \w{x)\^ '^■.nl+ UkUk+i ■■ 

k 

For the asymmetric parts, we first use Eq. (IA.9P obtaining 
gi dx : p^{x)pi{x) 



(A.15) 



gi_ 

8 



N 



dx \w{xr E 

and then finally the more complicated Eq. (jA.lOp 
dx : p^{xf ~ pi{xf 



2 \^'^fe+i^fe + '^k Cfe+1 



(A.16) 



52 

2 



92 

8 



N 

/ dx \w{x)\^ E : nfe(4cfc_i + 4-iCfc + c^^fe'+i + 4+iCfe) : (A-17) 



Introducing constants that parameterize the on-site interactions and the strength of 
the underlying lattice ([7]), our final Hamiltonian looks as follows 



Q k o 



: nkUk+i 



16 



Y H [("fe - l)4(cfe-i + Cfe+i) + H.C.] 



Y^icf^.cl+n.c.) 



(A.18) 



Completing terms and replacing the sum over k with a sum over nearest neighbors, 
we arrive at the desired model ([1]) with the parametrization given already in Eq. ((8|). 
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